Simulations of Interacting Membranes 
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The liquid crystalline model biomembrane system consisting of a stack of interacting membranes 
is studied by the newly developed Fourier Monte Carlo simulation technique. In comparison to 
perturbation theory, substantial quantitative discrepancies are found that affect determination of 
interbilayer interactions. A harmonic theory is also routinely used to interpret x-ray scattering 
line shapes; this is shown to be valid because the distance dependence of the simulated correlation 
functions can be fairly well fit by the harmonic theory. 

PACS numbers: 87.10. +e 02.70.Lq 61.30.Cz 
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Stacks of lipid bilayers (see Fig. are model systems 
for biomembranes that are much studied for two reasons. 
First, such stacks diffract fairly well and this facilitates 
determination of the structure of individual membranes, 
which is of primary interest in biophysics. However, these 
stacks are not crystals with the long range order that is 
assumed in traditional biophysical analysis, but smectic 
liquid crystals with quasi-long-range order. Therefore, 
quantitative use of the scattering intensity for structure 
determination requires correction for the fluctuations en- 
demic in such systems Q]. A harmonic theory that 
predicts power law tails in the scattering line shapes fits 
membrane data very well [^,^, but the anharmonicities 
that are inherent in realistic potentials have remained a 
concern for quantitative interpretation |6|j^ , even though 
a renormalization group analysis suggested that such ef- 
fects are small [0. 
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FIG. 1. Snapshot of a slice through a simulated stack of 
M = 8 two-dimensional L x L fluctuating membranes. Since 
internal membrane structure is irrelevant here, each mem- 
brane is depicted as a line. The average position of each 
membrane is shown by a dashed line. 



Stacks of bilayers are also much studied because they 
provide ideal environments to study fundamental inter- 
actions between bilayers, especially since the range of 
interbilayer distances a can be systematically varied by 
applying osmotic pressure P M . The corresponding the- 



ory 1 10 1 is an approximate first-order perturbation theory 
that again relies on harmonic assumptions, such as the 
normality of the probability distribution function (pdf) 
of the interbilayer spacing. While this theory has been 
a valuable guide, use of it to extract fundamental inter- 
bilayer interactions from P{a) data may be inaccurate. 

Both these issues are addressed using Monte Carlo sim- 
ulations with realistic intermembrane potentials for the 
biologically relevant regime where the interbilayer water 
spacing a is of order 5— 30 A and each membrane is flexible 
with bending modulus i^c~10 — 2bkBT . In this regime, 
sometimes called the soft confinement regime |Q, it is 
usually supposed that the primary interbilayer interac- 
tions for dipolar lipids are the attractive van der Waals 
potential and the repulsive hydration potential. 
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where z is the local distance between two membranes 
These interactions are significantly anharmonic, 
to the extent that the potential of a membrane mid- 
way between two neighboring membranes (at the dashed 
positions in Fig. |^) may have a maximum instead of 
a minimum. The contrasting regime, sometimes called 
the hard confinement regime, consists of only excluded 
volume or steric interactions between neighboring mem- 
branes. That regime is appropriate when a is of order 
lOOA because the hydration force is short range (A~2Jl). 
For hard confinement the effective interbilayer force is 
the entropic fluctuation pressure which decays as 
]l2t . Fluctuation forces also exist in the soft confinement 
regime, and are determined by our simulations. 

In addition to the interbilayer interactions in Eq.(^, 
the energy of each membrane includes a bending term, 
proportional to the square of the local curvature of the 
membrane. Let Um{x, y) be the local displacement of the 
TOth membrane from its average position as shown on 
Fig. Periodic boundary conditions are imposed in the 
plane of each membrane and also along the stack, so that 
um=uo- The membranes can collide, but cannot overlap, 
so that Um+i + a>Um, where a is the average distance 
between membranes. The Hamiltonian of the stack is 
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then: 
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m=0 
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{AUm)'^+V{Urn+l+a-Urn) 



dxdy (2) 



The simulation method, called the Fourier Monte 
Carlo method, was developed for single membranes be- 
tween hard walls and is easily extended to stacks of 
membranes. Each membrane in the stack is represented 
by a complex array of dimensions N x N oi Fourier dis- 
placement amplitudes. Instead of moving one lattice site 
at a time, moves are made in Fourier space and a whole 
membrane is displaced in each move. This allows larger 
moves and faster equilibration, without incurring large 
increases in the bending energy. One difference with our 
previous simulations M] is that a fixed osmotic pressure 
P ensemble is employed instead of the previous fixed a 
ensemble, so that a is obtained as a function of P rather 
than vice versa. Of course, use of the P ensemble is 
fundamentally no different, but it does have better con- 
vergence properties that we now discuss. 
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FIG. 2. Effect of finite density (TV = 6, ...,32) on 
a(solid squares) and A(open circles and right hand ordi- 
nate) for realistic interaction parameters given in and 
for P = W^erg/cm^. 



Simulations performed systematically as a function of 
lattice size, density of lattice points, and number of mem- 
branes in the stack show that accurate results for infinite, 
continuous membranes in infinite stacks can be obtained 
at one P in real time of the order of one day on a Pen- 
tium Pro PC. The most sensitive finite-size parameter is 
the "density" of each membrane N/L, since when N is 
varied from 6 to 32, the root mean square fluctuation in 
nearest neighbor spacing, defined as A, can easily change 
by 40%, and the changes in a are also significant. This is 
shown in Fig.^, which also shows that accurate values can 
be obtained by extrapolation. By comparison, variations 
with lateral system size L at fixed N/L are negligible 
(«0.2% for L>70QA), as are variations with M for 
M>8) for a and A. 

It may be noted that stacks of several membranes 
(Af«4) have been previously considered but 



mostly for the critical phenomenon of unbinding; this 
occurs in the limit of large average membrane spacing, 
where the van der Waals interaction is the main one, in 
addition to the spatial constraints. We have also per- 
formed simulations in the hard confinement regime and 
obtained results for the Helfrich fluctuation free energy 
CfiT"^ / KcO? with c/i«0.1, in agreement with [|l4|Jlf 
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FIG. 3. Simulation results (symbols) and perturbation the- 
ory (lines) for A and logP versus a for the parameter set in 
wA with attractive interaction (solid squares and lines) and 
with no attractive interaction (open circles and dashed lines). 



Figure ^ shows results for P and A as functions of 
a for realistic values of the potentials Also shown 

are results for a simpler case when the attractive van der 
Waals force is absent so that the potential experienced by 
each membrane has a minimum in the middle between its 
neighbors and is therefore more like a harmonic potential. 
Fig.^ also shows results based on the perturbation ap- 
proximation []lo|] which was developed for a single mem- 
brane between hard walls. After comparing to harmonic 
theory for multiple membranes we adjusted [|o| for 
the case of multiple membranes by putting a factor of 4/7r 
into the relations for the fluctuational free energy and for 
A^, and then followed essentially the same procedure as 
in [0. With this factor, agreement between the per- 
turbation theory and the simulations is quite good for 
small a where the total interaction is harmonic-like be- 
cause the repulsive hydration potential is dominant and 
the fluctuations are relatively small. The simulations and 
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the theory also agree for P{a) when there is no attrac- 
tive van der Waals interaction. However, at larger a, the 
simulated A increases with a faster than the perturba- 
tion approximation for either intermembrane potential. 
A large difference between the simulation results and the 
perturbation theory occurs when the potential has an 
attractive van der Waals part. The theory predicts that 
Oo = 25. SA for P = 0, more than 5A larger than the 
true value Qq = 20.2±0.1j4 obtained by simulation. It 
is also of interest to compare the values of /j, = (A/a)^ 
to the hard confinement values. In Fig.|| the range ofu 
is 0.06 — 0.12, in good agreement with experiment [gl, 
but considerably smaller than the hard confinement esti- 
mates 0.16 - 0.21 [RiillilMlil. 
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FIG. 4. Distribution of the nearest neighbor distance for 
parameters in and P = lO^erg/cm^ for different mem- 
brane densities. The bottom figure shows the same results as 
the top one, scaled and shifted so as to facilitate comparison 
to the standard normal distribution, denoted A'^(0, 1). 

Existing approximations assume that the dis- 

tribution of distance between two nearest neighbor mem- 
branes is essentially Gaussian. Simulation results in Fig. 
^ demonstrate that the actual pdf differs substantially 
from a normal distribution. Although the distribution 
is more non-normal for small N, it is nevertheless clear 
that even the limiting distribution is not normal. A dis- 
tinguishing feature of the true pdf is a rapid decay at 
small distances. This validates the use of a cutoff 
near 2 A to avoid the formal divergence of the vdW poten- 
tial. Another feature evident in Fig. || is the asymmetry 



of the actual pdf which shows that large fluctuations to 
larger intermembrane spacings are more probable than 
large fluctuations to smaller spacings. Overall, the shape 
of pdf is consistent with that of the interaction potential. 

We turn next to the issue of whether harmonic fluc- 
tuation theory should be expected to be reliable 
when interpreting detailed line shapes |^,|^ from stacks 
of bilayers and smectic liquid crystals that have strongly 
anharmonic interactions. The most important quantity 
for determining the line shapes for powder samples |^ is 
the correlation function in the z-direction, which is es- 
sentially the k dependence of the mean square relative 
displacement of two membranes of the stack A^(fc) = 
{{u{'m + fc) — u{m)Y), where m and m + k are the mem- 
brane indices, and averages over m are performed for 
simulation efficiency. Fig.H shows profiles of A(fc), ob- 
tained for stacks with various numbers M of membranes. 
Convergence with M suggests that values of A(fc) are 
sufficiently accurate for k < M /A. However, to minimize 
the finite size effect in a comparison with harmonic the- 
ory. Fig. |5| compares the results of the simulation with 
M = 32 with the exact harmonic result, also for M = 32. 
In the harmonic theory the bare interbilayer interactions 
are approximated with a compression modulus B. In Fig. 
H a value of i? = 1.9-\Q^^erg/crn^ was chosen to match 
the large k end of the M = 32 curve. The resulting A(fc) 
profile allows one to see that A=A(1) is in fact a good 
proxy for describing the long-range correlations, since the 
difference between A, implied by the "harmonic curve" , 
and the actual A for the stack is about only 0.2A, i.e. 
relatively small compared to A = 4.6A. Another way to 
see how interactions are renormalized from short to long 
range is to compute, for different fc, the implied strength 
B{k) of the harmonic potential that would result in the 
same value of A(fc) as obtained by simulations for a stack 
with realistic interactions. The bottom panel of Fig.^ 
presents a plot of the harmonic value of B{k) required to 
give the simulation value of A(/c). This shows that, for 
large k, the system can be reasonably well approximated 
by one with harmonic interactions with constant B. 

How is it that the harmonic theory works quite well 
for the correlation functions in the preceding paragraph 
and not so well for P and A in Fig. |[? The answer 
is that the perturbation theory does not yield the best 
value of B\ for the example in Fig.^ the theory yields 
a larger B = 5A-10^^erg/cm^ which accounts for the 
smaller value of A in Fig. 

We turn flnally to the entropic fluctuation pressure in 
a stack of membranes, which is defined to be the differ- 
ence between the applied pressure, and the pressure due 
to the bare van der Waals and hydration interactions. 
Perturbation theory |]l0| , experiment |^ and simulations 
on a single membrane between hard walls ||l^] all agree 
that the decay of the fluctuation pressure is closer to ex- 
ponential, with a decay length A , although the value of 
Xfi found in both experiment and simulations is larger 
than the perturbation theory prediction Xfi — 2A. Fig. 
^ shows the same result for simulations of stacks. 
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FIG. 5. Root mean square fluctuations A(fc) between fcth 
neighbor membranes for a stack with difi'erent numbers M 
of membranes and for the parameter set in jl^ (except that 
L = 1400A) at P = 3.16- 10'^ erg/cm^. Also shown is A(fc), 
exactly computed for the case of harmonic interactions with 
compression modulus B = 1.9- 10^"^ erg/cm*. The bottom fig- 
ure shows the effective harmonic compression modulus B{k) 
for each k and M. 
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FIG. 6. Simulation results for Pfi vs. a for the parameter 
set in and also for H = 0. The slope is A/; = 4.34A. 



A long range goal is to obtain values of the interbilayer 
interaction parameters; the traditional analysis |^ uses 
osmotic pressure P{a) data, which has recently been sup- 
plemented by fluctuation A(a) data ||^. One of the main 
results of this paper indicates that the A data are indeed 



valid, even though the analysis of the basic x-ray scat- 
tering data is based on a harmonic theory. However, the 
intrinsic anharmonic nature of realistic interactions be- 
tween bilayers in stacks makes it difficult to devise quan- 
titatively accurate analytic or perturbation theories. We 
show here that the Fourier Monte Carlo method is sufR- 
ciently fast that it provides a viable alternative. Indeed, 
it is now possible to consider using it as part of a com- 
prehensive data analysis program to determine the best 
values of the fundamental interaction parameters. 
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